/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */

/*
 * This sample implements Mersenne Twister random number generator 
 * and Cartesian Box-Muller transformation on the GPU.
 * See supplied whitepaper for more explanations.
 */


// Utilities and system includes
//#include "shrUtils.h"
//#include "shrQATest.h"
//#include <cutil_inline.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "MersenneTwister.h"

///////////////////////////////////////////////////////////////////////////////
// Common host and device function 
///////////////////////////////////////////////////////////////////////////////
//ceil(a / b)
extern "C" int iDivUp(int a, int b){
    return ((a % b) != 0) ? (a / b + 1) : (a / b);
}

//floor(a / b)
extern "C" int iDivDown(int a, int b){
    return a / b;
}

//Align a to nearest higher multiple of b
extern "C" int iAlignUp(int a, int b){
    return ((a % b) != 0) ?  (a - a % b + b) : a;
}

//Align a to nearest lower multiple of b
extern "C" int iAlignDown(int a, int b){
    return a - a % b;
}

///////////////////////////////////////////////////////////////////////////////
// Reference MT front-end and Box-Muller transform
///////////////////////////////////////////////////////////////////////////////
//extern "C" void initMTRef(const char *fname);
//extern "C" void RandomRef(float *h_Random, int NPerRng, unsigned int seed);
//extern "C" void BoxMullerRef(float *h_Random, int NPerRng);

///////////////////////////////////////////////////////////////////////////////
// Fast GPU random number generator and Box-Muller transform
///////////////////////////////////////////////////////////////////////////////
#include "MersenneTwister_kernel.cu"

///////////////////////////////////////////////////////////////////////////////
// Data configuration
///////////////////////////////////////////////////////////////////////////////
#define     SUM_N 4096
const int    PATH_N = 96000000; //Quantidade de Números a ser gerados.
int N_PER_RNG = iAlignUp(iDivUp(PATH_N, MT_RNG_COUNT), 2); //Quantidade de Números por RING
int    RAND_N = MT_RNG_COUNT * N_PER_RNG;

unsigned int SEED = 777;

__global__ void sum_vector(float *A, float *B, float *C)
{
	int id = blockDim.x * blockIdx.x + threadIdx.x;

	if(id < SUM_N){
		C[id] = A[id] + B[id];
	}
}

__global__ void sum_vector_stream(float *A, float *B, float *C)
{
	int id = blockDim.x * blockIdx.x + threadIdx.x, t;

	for(t=0;t<1024;t++){

		if(id < SUM_N){
			C[id] = A[id] + B[id +(t*MT_RNG_COUNT)];
		}

	}
}

//#define DO_BOXMULLER
///////////////////////////////////////////////////////////////////////////////
// Main program
///////////////////////////////////////////////////////////////////////////////
int main(int argc, char **argv)
{
    // Start logs
    //shrQAStart(argc, argv);

    //shrSetLogFileName ("MersenneTwister.txt");
    printf("%s Starting...\n\n", argv[0]);

    float *d_Rand, *Vector_C, *Vector_A, *Vector_B,
        /*  *h_RandCPU,*/
          *h_RandGPU, *vector_GPU, *vector, *HostRandGPU, *DevPointer;

    double /*rCPU, rGPU, delta, sum_delta, max_delta, sum_ref, L1norm,*/ gpuTime;
    
    cudaEvent_t start, stop;

    int i, j, ver = 1;
    float hTimer;
    
    cudaDeviceProp prop;
    int whichDevice;
    
    cudaGetDevice( &whichDevice );
    cudaGetDeviceProperties( &prop, whichDevice );
    
    if( prop.canMapHostMemory != 1){
    
    	printf("Device cannot map memory.\n");
    	return 0;
    
    }
    if( !prop.deviceOverlap ){
    
    	printf("Device will not handle overlaps, so no speed up streams.\n");
    	return 0;
    
    }
    
    printf("Device - %s\n\n", prop.name);
    
    //if( cutCheckCmdLineFlag(argc, (const char**)argv, "device") )
    //    cutilDeviceInit(argc, argv);
    //else
    //cudaSetDevice( cutGetMaxGflopsDeviceId() );

    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    printf("Initializing data for %i samples...\n", PATH_N);
        //h_RandCPU  = (float *)malloc(RAND_N * sizeof(float));
	vector_GPU = (float *)malloc(SUM_N * sizeof(float));
	vector = (float *)malloc(SUM_N * sizeof(float));
        //h_RandGPU  = (float *)malloc(RAND_N * sizeof(float));
        cudaHostAlloc((void **)&HostRandGPU, RAND_N * sizeof(float), cudaHostAllocWriteCombined | cudaHostAllocMapped );
        cudaHostGetDevicePointer( &DevPointer, HostRandGPU, 0);
        cudaMalloc((void **)&d_Rand, RAND_N * sizeof(float));
        cudaMalloc((void **)&Vector_C, SUM_N * sizeof(float));
	cudaMalloc((void **)&Vector_A, SUM_N * sizeof(float));
	cudaMalloc((void **)&Vector_B, SUM_N * sizeof(float));
	
	for(i=0; i<SUM_N;i++){
		vector[i] = 7;
	}

	cudaMemcpy(Vector_A, vector, SUM_N * sizeof(float), cudaMemcpyHostToDevice);

    printf("Loading GPU twisters configurations...\n");
        //const char *raw_path = shrFindFilePath("MersenneTwister.raw", argv[0]);
        const char *dat_path = "/home/emeson/NVIDIA_CUDA-5.0_Samples/EmesonSantana/MersenneTwister/data/MersenneTwister.dat";//shrFindFilePath("MersenneTwister.dat", argv[0]);
        //initMTRef(raw_path);
        loadMTGPU(dat_path);
        seedMTGPU(SEED);

	//int numIterations = 10;
	//for (i = -1; i < numIterations; i++)
	//{
		//if (i == 0) 
		//{

/////////////////////////////////////////////////////////////////////////////////////////////////////////
// Mersenne Twister GPU Test 1*
/////////////////////////////////////////////////////////////////////////////////////////////////////////
	printf("Generating random numbers on GPU and Summing the vectors ***Test 1***...\n\n");
	//cutilSafeCall( cutilDeviceSynchronize() );
	cudaSetDeviceFlags( cudaDeviceMapHost );
	
	cudaEventRecord( start, 0);
		//}
        RandomGPU<<<32, 128>>>(DevPointer, N_PER_RNG); //
        //cutilCheckMsg("RandomGPU() execution failed\n");
//#ifdef DO_BOXMULLER        
        //BoxMullerGPU<<<32, 128>>>(d_Rand, N_PER_RNG);
        //cutilCheckMsg("BoxMullerGPU() execution failed\n");
//#endif
    //}
    //cutilSafeCall( cutilDeviceSynchronize() );

	//shrLog("\nReading back the results in Host and sum vectors...\n");
        	//cudaMemcpy(h_RandGPU, d_Rand, RAND_N * sizeof(float), cudaMemcpyDeviceToHost); // Números aleatórios Device To Host
	for(j=0; j< RAND_N/SUM_N;j++){
		//cudaMemcpy(Vector_B, &h_RandGPU[j*SUM_N], SUM_N * sizeof(float), cudaMemcpyHostToDevice);//Números aleatórios Host To Device de SUM_N de cada vez
		sum_vector<<<32, 128>>>(Vector_A, &DevPointer[j*SUM_N], Vector_C);// SOMA dos Vetores de SUM_N de cada vez paralelamente
		//cutilSafeCall( cutilDeviceSynchronize() );
	}
	
	cudaThreadSynchronize();

    cudaEventRecord( stop, 0);
    cudaEventSynchronize( stop );
    
    cudaEventElapsedTime( &hTimer, start, stop);

    gpuTime = 1.0e-3 * hTimer;// /(double)numIterations;
    printf("MersenneTwister 1*, Time = %.5f s, Size = %u Numbers Generated - %u Numbers Used - %u Numbers Per Time in Parallel Sum, Requisitions Numbers = %u Times, NumDevsUsed = %u, Numbers of Blocks in MersenneGPU = %u, Numbers of Threds Per Block in MersenneGPU = %u, Numbers of Blocks in Parallel Sum = %u, Numbers of Threds Per Block in Parallel Sum = %u\n\n", 
               gpuTime, RAND_N, RAND_N, SUM_N, RAND_N/SUM_N, 1, 32, 128, 32, 128);

/////////////////////////////////////////////////////////////////////////////////////////////////////////
// Mersenne Twister GPU Test 2*
/////////////////////////////////////////////////////////////////////////////////////////////////////////
	printf("Generating random numbers on GPU and Summing the vectors ***Test 2***...\n\n");
	//cutilSafeCall( cutilDeviceSynchronize() );
	cudaEventRecord( start, 0);
		//}
        RandomGPU<<<32, 128>>>(d_Rand, N_PER_RNG); //
        //cutilCheckMsg("RandomGPU() execution failed\n");
//#ifdef DO_BOXMULLER        
        //BoxMullerGPU<<<32, 128>>>(d_Rand, N_PER_RNG);
        //cutilCheckMsg("BoxMullerGPU() execution failed\n");
//#endif
    //}
    //cutilSafeCall( cutilDeviceSynchronize() );

	//shrLog("\nReading back the results in Host and sum vectors...\n");
	for(j=0; j< RAND_N/SUM_N;j++){
		sum_vector<<<32, 128>>>(Vector_A, &d_Rand[j*SUM_N], Vector_C);// SOMA dos Vetores de SUM_N de cada vez paralelamente
		//cutilSafeCall( cutilDeviceSynchronize() );
	}

    cudaEventRecord( stop, 0);
    cudaEventSynchronize( stop );
    
    cudaEventElapsedTime( &hTimer, start, stop);

    gpuTime = 1.0e-3 * hTimer;// /(double)numIterations;
    printf("MersenneTwister 2*, Time = %.5f s, Size = %u Numbers Generated - %u Numbers Used - %u Numbers Per Time in Parallel Sum, Requisitions Numbers = %u Times, NumDevsUsed = %u, Numbers of Blocks in MersenneGPU = %u, Numbers of Threds Per Block in MersenneGPU = %u, Numbers of Blocks in Parallel Sum = %u, Numbers of Threds Per Block in Parallel Sum = %u\n\n", 
               gpuTime, RAND_N, RAND_N, SUM_N, RAND_N/SUM_N, 1, 32, 128, 32, 128);


/////////////////////////////////////////////////////////////////////////////////////////////////////////
// Mersenne Twister GPU Test 3*
///////////////////////////////////////////////////////////////////////////////////////////////////////// 
	printf("Generating random numbers on GPU and Summing the vectors ***Test 3***...\n\n");
	//cutilSafeCall( cutilDeviceSynchronize() );
	N_PER_RNG = 1;
	cudaFree(d_Rand);
	h_RandGPU  = (float *)malloc(MT_RNG_COUNT * sizeof(float));
        cudaMalloc((void **)&d_Rand, MT_RNG_COUNT * sizeof(float));

	cudaEventRecord( start, 0);

        // Gerando Números aleatórios
        //cutilCheckMsg("RandomGPU() execution failed\n");
    	//cutilSafeCall( cutilDeviceSynchronize() );

	//shrLog("\nReading back the results in Device and sum vectors...\n");
	for(j=0; j< RAND_N/SUM_N;j++){
		RandomGPU<<<32, 128>>>(d_Rand, N_PER_RNG);
		sum_vector<<<32, 128>>>(Vector_A, d_Rand, Vector_C); // SOMA dos Vetores de SUM_N de cada vez paralelamente
		//cutilSafeCall( cutilDeviceSynchronize() );
	}

    cudaEventRecord( stop, 0);
    cudaEventSynchronize( stop );
    
    cudaEventElapsedTime( &hTimer, start, stop);

	cudaMemcpy(vector_GPU, Vector_C, SUM_N * sizeof(float), cudaMemcpyDeviceToHost);
    gpuTime = 1.0e-3 * hTimer;// /(double)numIterations;
    printf("MersenneTwister 3*, Time = %.5f s, Size = %u Numbers Generated - %u Numbers Used - %u Numbers Per Time in Parallel Sum, Requisitions Numbers = %u Times, NumDevsUsed = %u, Numbers of Blocks in MersenneGPU = %u, Numbers of Threds Per Block in MersenneGPU = %u, Numbers of Blocks in Parallel Sum = %u, Numbers of Threds Per Block in Parallel Sum = %u\n\n", 
               gpuTime, RAND_N, RAND_N, SUM_N, RAND_N/SUM_N, 1, 32, 128, 32, 128);

/////////////////////////////////////////////////////////////////////////////////////////////////////////
// Mersenne Twister GPU Test 4*
/////////////////////////////////////////////////////////////////////////////////////////////////////////

	printf("Generating random numbers on GPU and Summing the vectors ***Test 4***...\n\n");

	cudaStream_t stream_mt, stream_sum;
	cudaStreamCreate( &stream_mt );
	cudaStreamCreate( &stream_sum );
	
	N_PER_RNG = 1024;

	cudaFree(d_Rand);
	cudaMalloc((void **)&d_Rand, 2 * 1024 * MT_RNG_COUNT * sizeof(float));

	cudaEventRecord( start, 0);
	
	for(j=0; j< (RAND_N/(SUM_N*1024)) +1;j++){
		
		RandomGPU<<<32, 128, 0, stream_mt>>>(&d_Rand[(j%2)*1024 * MT_RNG_COUNT], N_PER_RNG);
		cudaStreamSynchronize( stream_mt );
		sum_vector_stream<<<32, 128, 0, stream_sum>>>(Vector_A, &d_Rand[(j%2)*1024 * MT_RNG_COUNT], Vector_C); // SOMA dos Vetores de SUM_N de cada vez paralelamente
	}
	
	cudaStreamSynchronize( stream_sum );
	cudaStreamSynchronize( stream_mt );

    	cudaEventRecord( stop, 0);
    	cudaEventSynchronize( stop );
    
    	cudaEventElapsedTime( &hTimer, start, stop);

	cudaMemcpy(vector_GPU, Vector_C, SUM_N * sizeof(float), cudaMemcpyDeviceToHost);
    gpuTime = 1.0e-3 * hTimer;// /(double)numIterations;
    printf("MersenneTwister 4*, Time = %.5f s, Size = %u Numbers Generated - %u Numbers Used - %u Numbers Per Time in Parallel Sum, Requisitions Numbers = %u Times, NumDevsUsed = %u, Numbers of Blocks in MersenneGPU = %u, Numbers of Threds Per Block in MersenneGPU = %u, Numbers of Blocks in Parallel Sum = %u, Numbers of Threds Per Block in Parallel Sum = %u\n\n", 
               gpuTime, RAND_N, (j+1)*SUM_N*1024, SUM_N, RAND_N/SUM_N, 1, 32, 128, 32, 128);


////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Checking
////////////////////////////////////////////////////////////////////////////////////////////////////////////
    printf("Checking GPU results...\n");
	cudaMemcpy(vector_GPU, Vector_C, SUM_N * sizeof(float), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_RandGPU, &d_Rand[((j%2)* 1024 * MT_RNG_COUNT) + (1023 * MT_RNG_COUNT)], MT_RNG_COUNT * sizeof(float), cudaMemcpyDeviceToHost);
	for(i= 0; i < SUM_N;i++){
		if(vector_GPU[i] - h_RandGPU[i] < vector[i] - 1.0e-6) ver = 0;
	}
        //shrLog(" ...generating random numbers on CPU using reference generator\n");
        //RandomRef(h_RandCPU, N_PER_RNG, SEED);

        //#ifdef DO_BOXMULLER
            //shrLog(" ...applying Box-Muller transformation on CPU\n");
            //BoxMullerRef(h_RandCPU, N_PER_RNG);
        //#endif

      /*shrLog(" ...comparing the results\n\n");
        max_delta = 0;
        sum_delta = 0;
        sum_ref   = 0;
        for(i = 0; i < MT_RNG_COUNT; i++)
            for(j = 0; j < N_PER_RNG; j++){
                rCPU = h_RandCPU[i * N_PER_RNG + j];
                rGPU = h_RandGPU[i + j * MT_RNG_COUNT];
                delta = fabs(rCPU - rGPU);
                sum_delta += delta;
                sum_ref   += fabs(rCPU);
                if(delta >= max_delta) max_delta = delta;
            }
    L1norm = (float)(sum_delta / sum_ref);
    shrLog("Max absolute error: %E\n", max_delta);
    shrLog("L1 norm: %E\n\n", L1norm);*/

    printf("Shutting down...\n");
        cudaFree(d_Rand);
	cudaFree(Vector_A);
	cudaFree(Vector_B);
	cudaFree(Vector_C);
	cudaFreeHost(HostRandGPU);
        free(h_RandGPU);
	free(vector_GPU);
	free(vector);
	cudaEventDestroy( start );
	cudaEventDestroy( stop );
	cudaStreamDestroy( stream_mt );
	cudaStreamDestroy( stream_sum );
        //free(h_RandCPU);
	
   ver ?/*(L1norm < 1e-6) ?*/ printf("PASSED\n") : printf("FAILED\n");
}
